The TESS-Keck Survey. XII. A Dense 1.8 R_⊕ Ultra-short-period Planet Possibly Clinging to a High-mean-molecular-weight Atmosphere after the First Gigayear

TOI-1347 (TIC 229747848) is a late G-type (T_eff = 5464 ± 100 K) star that is relatively active ($\mathrm{log}{R}_{{\rm{HK}}}^{{\prime} }=-4.66$), showing strong variability in both photometry (see Section 4) and RVs (see Section 5).

We obtained high-resolution spectra of TOI-1347 with HIRES (with the B3 decker, R = 67,000) and the High Accuracy Radial Velocity Planet Searcher-North (HARPS-N; R = 115,000); see Section 5 for details. We applied the SpecMatch-Synthetic (Petigura 2015) algorithm to the HIRES spectrum to measure T_eff, M_⋆, Fe/H, and $\mathrm{log}g$. These parameters, as well as the J and K magnitudes from the TICv8 catalog (Paegert et al. 2021) and the Gaia parallax (Gaia Collaboration et al. 2023), were input into isoclassify (Huber et al. 2017; Berger et al. 2020) to constrain R_⋆. We also used KeckSpec (A. Polanski 2024, in preparation) to compute alpha elemental abundances. Stellar parameters for TOI-1347 are included in the catalog of MacDougall et al. (2023). We verified that our values agree with the values therein to within 1σ.

We similarly derived spectroscopic parameters from the coadded HARPS-N spectra of TOI-1347 using the FASMA spectral synthesis package (Tsantaki et al. 2018, 2020). We verified that all derived quantities were in agreement with the HIRES values. In particular, the HARPS-N data are of sufficient resolution to detect $v\sin {i}_{\star }$, whereas with HIRES we only obtain an upper limit. Table 1 lists our adopted stellar parameters.

Table 1. Stellar Parameters of TOI-1347

Parameter	Value	Unit	Source
TIC ID	TIC 229747848	⋯	TIC v8.2 a
R.A.	18:41:18.4	hh:mm:ss	TIC v8.2 a
Decl.	+70:17:24.19	dd:mm:ss	TIC v8.2 a
V magnitude	11.168 ± 0.013	⋯	TIC v8.2 a
TESS magnitude	10.7157 ± 0.0066	⋯	TIC v8.2 a
J magnitude	10.011 ± 0.02	⋯	TIC v8.2 a
K magnitude	9.616 ± 0.015	⋯	TIC v8.2 a
Gaia magnitude	11.2076 ± 0.0005	⋯	Gaia DR3 b
Parallax	6.7803 ± 0.0112	mas	Gaia DR3 b
R.A. proper motion	−6.0883 ± 0.0150	mas yr−1	Gaia DR3 b
Decl. proper motion	26.1020 ± 0.0150	mas yr−1	Gaia DR3 b
Galactic U	−5.31 ± 0.09	km s−1	Derived from Gaia DR3 astrometry c
Galactic V	−7.72 ± 0.46	km s−1	Derived from Gaia DR3 astrometry c
Galactic W	4.91 ± 0.23	km s−1	Derived from Gaia DR3 astrometry c
Luminosity	0.55 ± 0.02	L⊙	Isoclassify d
Radius	0.83 ± 0.03	R⊙	Isoclassify d
Mass	0.913 ± 0.033	M⊙	SpecMatch-Synthetic e
Teff	5464 ± 100	K	SpecMatch-Synthetic e
$\mathrm{log}g$	4.64 ± 0.10	⋯	SpecMatch-Synthetic e
[Fe/H]	0.04 ± 0.06	dex	SpecMatch-Synthetic e
[α/Fe]	−0.03 ± 0.06	dex	KeckSpec f
$v\sin {i}_{\star }$	<3	km s−1	SpecMatch-Synthetic e
$v\sin {i}_{\star }$	2.9 ± 0.1	km s−1	FASMA g
veq	2.61 ± 0.11	km s−1	From Prot and R*
$\mathrm{log}{R}_{{\rm{HK}}}^{{\prime} }$	−4.66 ± 0.05	⋯	From Keck/HIRES spectrum
Prot	16.1 ± 0.3	days	From ACF of TESS photometry (Section 4.1)
Prot	16.3 ± 0.6	days	From TESS-SIP with TESS photometry (Section 4.1)
Prot	16.2 ± 0.3	days	From GP fit to RVs (Section 5.1)
Age	1.7 ± 0.1	Gyr	From Prot using empirical relations of Bouma et al. (2023)
Age	1.33 ± 0.06	Gyr	From Prot using empirical relations of Mamajek & Hillenbrand (2008)
Age	1.6 ± 0.4	Gyr	From $\mathrm{log}{R}_{{\rm{HK}}}^{{\prime} }$ using empirical relations of Mamajek & Hillenbrand (2008)
Age	${0.8}_{-0.6}^{+1.1}$	Gyr	From isochronal fitting by MacDougall et al. (2023)
Age	>650	Myr	From Li and comparison to Hyades (Berger et al. 2018)
Age	1.4 ± 0.4	Gyr	Adopted value

Notes.

^a Paegert et al. (2021). ^b Gaia Collaboration et al. (2023). ^c Using local standard of rest from Coşkunoǧlu et al. (2011), as implemented in PyAstronomy.pyasl.gal_uvw (Czesla et al. 2019). ^d Huber et al. (2017), using the SpecMatch-Synthetic results from an HIRES spectrum. Uncertainties are random errors for the adopted model grid. To account for systematic errors, see Tayar et al. (2022). ^e Petigura (2015), using an HIRES spectrum taken at R ∼ 67,000 with no iodine. ^f A. Polanski (2024, in preparation), using a HIRES spectrum taken at R ∼ 67,000 with no iodine. ^g Tsantaki et al. (2018, 2020), using the coadded HARPS-N spectra.

Download table as:  ASCIITypeset image

We constrained the age of the host star using several methods. Using a stellar rotation period of 16.1 ± 0.3 days (see Section 4.1 for details), we estimated the gyrochronological age of the star. The scaling relation of Mamajek & Hillenbrand (2008) gives an age of 1.33 ± 0.06 Gyr. Using the latest empirical relations of Bouma et al. (2023), the age is 1.7 ± 0.1 Gyr. We also obtained an age estimate using chromospheric activity in the Ca ii H and K lines from our HIRES spectra. We measured $\mathrm{log}{R}_{\mathrm{HK}}^{{\prime} }=-4.66\pm 0.05$ using the method of Isaacson & Fischer (2010). Combining with the calibration of Mamajek & Hillenbrand (2008), the corresponding age is 1.6 ± 0.4 Gyr. We also looked for the lithium doublet in our HIRES spectra but were unable to detect the lithium doublet above the noise floor in the continuum (see Figure 1). We placed an upper limit of 2 m Å (95% confidence) on the equivalent width. According to Berger et al. (2018), the star is consistent with field stars and is most likely older than the Hyades cluster (∼650 Myr). Lastly, MacDougall et al. (2023) derived an age of ${0.8}_{-0.6}^{+1.1}$ Gyr from isochronal fitting.

Zoom In Zoom Out Reset image size

The age uncertainties above do not account for systematic errors. Therefore, we combined age indicators with an unweighted mean. We adopted a wide age uncertainty of 0.4 Gyr to reflect the systematic uncertainties between the different age estimators. Our best age estimate for TOI-1347 is thus 1.4 ± 0.4 Gyr. Importantly, this age is longer than the timescale on which photoevaporation operates, which is typically confined to the first few hundred megayears when the star is active in X-rays and extreme UV (Ribas et al. 2005; Tu et al. 2015).

The effect of a starspot on photometry is to reduce the observed (integrated) intensity when on the visible hemisphere. The net effect of many spots is quasiperiodic variability that can be treated as time-correlated noise (Haywood et al. 2014; Rajpaul et al. 2015; Aigrain & Foreman-Mackey 2023). The TESS photometry of TOI-1347 shows strong rotational variability with maxima/minima occurring roughly every ∼8 days (see the top panel in Figure 3). A Lomb–Scargle periodogram of the photometry shows a strong peak at 8 days (Figure 4). However, a closer examination of the light curve reveals that the depths of adjacent maxima/minima are dissimilar. In fact, the depths of alternating maxima/minima tend to have similar amplitudes. This can be explained if the star has a 16 day rotation period and multiple spot groups concentrated on opposite hemispheres of the star, an effect noted for a number of other stars (Holcomb et al. 2022). In fact, a periodogram analysis of the RV data set shows a strong peak at 16 days (Figure 4). Spots affect RVs in much the same way as photometry by breaking the flux balance across star's rotational velocity profile (Saar & Donahue 1997), in addition to suppressing the convective blueshift (Haywood et al. 2016).

Zoom In Zoom Out Reset image size

While the periodogram is essentially a Fourier decomposition showing the amplitude of the best-fitting sine wave at all possible periods, the autocorrelation function (ACF) measures the self-similarity of a time series as a function of time delay. Thus, if the variability in a time series is not sinusoidal (asymmetric, more complex shape), then the ACF will give a better estimate of the periodicity than the periodogram. The ACF of the TESS photometry has its highest peak at ∼16 days (Figure 3) with adjacent ACF peaks alternating between high and low amplitude. An analysis of the ACF using spinspotter (Holcomb et al. 2022) successfully identified the half-period effect by checking that the odd peaks are less than 10% the height of the even peaks. spinspotter returns a rotation period of P_rot = 16.1 ± 0.3 days by averaging the locations of the even-numbered peaks only. We also measured the stellar rotation period using the TESS Systematics Insensitive Periodogram (TESS-SIP) algorithm (Hedges et al. 2020). TESS-SIP simultaneously corrects for instrument systematics while performing the periodogram search on the TESS Simple Aperture Photometry. Using TESS-SIP for sectors 14–26 for TOI-1347, we measured a stellar rotation period of 16.34 ± 0.57 days.

With all of this considered, we adopt the ∼16 day solution as the rotation period of TOI-1347 and explain the Lomb–Scargle periodogram peak at P_rot/2 as arising from antipodal spot groups. The persistence of repeated peaks every 8 days in the ACF, even out to beyond 100 days, can be explained if the starspots on TOI-1347 live for many rotation periods. Following the prescription of Giles et al. (2017), we fit the observed ACF with an underdamped simple harmonic oscillator, including a second component with power at P/2:

$$\begin{eqnarray}&&y={e}^{-{\rm{\Delta }}t/\tau }\left[A\cos \left(\displaystyle \frac{2\pi {\rm{\Delta }}t}{P}\right)+B\cos \left(\displaystyle \frac{2\pi {\rm{\Delta }}t}{P/2}\right)+{y}_{0},\right].\end{eqnarray} \tag{ 1 }$$

In Equation (1), y is the ACF strength, and the coefficients A and B give the relative strengths of the antipodal spot groups. We used scipy.optimize (Virtanen et al. 2020) to find the maximum a posteriori (MAP) solution, then used that as a seed for a Markov Chain Monte Carlo (MCMC) analysis using emcee (Foreman-Mackey et al. 2013). The fit recovered P = 16.0 days, and the best-fit exponential decay timescale was τ = 45 days, roughly 3 times the rotation period. Both coefficients A = 0.05 and B = 0.18 were constrained to nonzero values, again supporting the multiple spot group hypothesis.

Lastly, our SpecMatch-Synthetic analysis of the HIRES spectrum in Section 2.1 yielded only an upper limit corresponding to the line spread function of HIRES (roughly 2.2 km s⁻¹). Masuda et al. (2022) recently showed that on the population level, such nondetections of $v\sin {i}_{\star }$ are most consistent with a <3 km s⁻¹ upper limit. Combining P_rot = 16.1 days with the stellar radius from Section 2.1, we get an equatorial rotational velocity of v_eq = 2.6 ± 0.1 km s⁻¹. This is consistent with a <3 km s⁻¹ upper bound for HIRES though it is slightly smaller than the measured HARPS-N value. Were the rotation period 8 days, we would instead have v_eq ∼ 5 km s⁻¹, which (barring a misaligned stellar inclination of ≲30°) would be detectable in the HIRES spectrum and inconsistent with the HARPS-N measurement. If astrophysical, the slightly larger $v\sin {i}_{\star }$ could be explained by differential rotation and long-lived starspots at higher latitudes.

We searched the TESS light curve for planetary transits with a Box-fitting Least Squares (BLS; Kovács et al. 2002) algorithm as described in Dai et al. (2021). We recovered the two planet candidates reported by ExoFOP at 0.84 and 4.84 days. We did not find any other transit signals with signal-to-noise ratio (S/N) > 6.5.

Our transit analysis closely follows that described in Dai et al. (2021). Briefly, we generated transit signals using the Python package Batman (Kreidberg 2015), parameterized with the stellar density (ρ_*) to break the degeneracy between the scaled semimajor axes (a/R_*) and impact parameters (b) of the transiting planets (Seager & Mallén-Ornelas 2003). We imposed a prior using the best-fit stellar density from our adopted M_* and R_* of ρ_* = 2.25 ± 0.26 g cm⁻³. For the limb darkening coefficients, we used the prior and parameterization scheme of Kipping (2013) for a quadratic limb darkening law (q₁ and q₂). The other transit parameters included the orbital period (P_orb), time of conjunction (T_c), planet-to-star radius ratios (R_p/R_⋆), scaled orbital distances (a/R_⋆, computed from stellar density and orbital period), orbital inclinations (cosi), orbital eccentricities (e), and the arguments of pericenter (ω). We initially allowed nonzero eccentricities for both TOI-1347 b and c; however, the posterior distributions are fully consistent with circular orbits. Neither the existing transit nor RV data (Section 5) support the detection of nonzero eccentricities. In our final fits, we chose to restrict both planets to circular orbits to reduce model complexity.

Our transit fitting pipeline takes the following steps:

• 1.  
  Fit a global model of all transit epochs, assuming no transit timing variations (TTVs). This model is found by maximizing the likelihood with the Levenberg-Marquardt method implemented in Python package lmfit (Newville et al. 2014).
• 2.  
  Search for TTVs. We held the transit shape parameters fixed and fit for the mid-transit times of each transit epoch. We removed out-of-transit variations with a quadratic model. We did not detect significant TTVs for either planet, so we continued with fixed orbital periods.
• 3.  
  Full model. We sampled the posterior distributions of both planets jointly using the MCMC framework implemented in emcee (Foreman-Mackey et al. 2013) with 128 walkers initialized near the MAP solution. We ran emcee for 50,000 steps, ensuring that this was many times longer than the autocorrelation of various parameters (typically 100 s of steps).

Figure 5 shows the phase-folded and binned transits of TOI-1347 b and c with the MAP model. Using the stellar radius derived in Section 2, we derived R_p,b = 1.8 ± 0.1 R_⊕ and R_p,c = 1.6 ± 0.1 R_⊕, in agreement with the radii measured by MacDougall et al. (2023;${R}_{p,b}={1.81}_{-0.06}^{+0.09}\,{R}_{\oplus }$ and ${R}_{p,c}\,={1.68}_{-0.07}^{+0.09}$).

Zoom In Zoom Out Reset image size

We searched the TESS light curve for phase-curve variations and secondary eclipses of TOI-1347 b. First, we masked the in-transit data and removed long-term variability (stellar and/or instrumental) using the method of Sanchis-Ojeda et al. (2013). This involves fitting a linear function of time to the out-of-transit data points within a window of 2× the orbital period, then dividing the best-fit function within that window, repeating for every data point in the light curve. This detrended light curve is then phase folded to the orbital period of TOI-1347 b.

The resultant phase curve is shown in Figure 6. We were able to detect a tentative phase-curve variation (3σ) and a secondary eclipse (2σ) using a joint model. To model the secondary eclipse, we simply modified the best-fit transit model by shifting the mid-transit times by half the orbital period (i.e., assuming e = 0) and turning off limb darkening (q₁ = q₂ = 0). The secondary eclipse depth (${\delta }_{\sec }$) is allowed to vary freely to account for a combined effect of reflected stellar light and thermal emission from the planet's nightside. For the phase curve, we used a Lambertian disk model (see, e.g., Demory et al. 2016) parameterized by an amplitude A and phase offset of the peak θ. We sampled the posterior distribution with an MCMC analysis similar to that described in Section 4.3. We found a secondary eclipse depth of ${\delta }_{\sec }=26\pm 12\,\mathrm{ppm}$ and a phase-curve amplitude of A = 28 ± 9 ppm. The peak of the phase-curve variation is shifted by 33° ± 14° to the west of the planet.

Zoom In Zoom Out Reset image size

The lower panel of Figure 6 compares the amount of thermal emission versus reflected light as a function of the planet's Bond albedo (A_B ). At high albedo, the planet is more reflective, and the equilibrium temperature will be lower. In this case, the phase curve will be dominated by reflected stellar light rather than thermal emission. This is necessary to explain the large secondary eclipse depth measured in the TESS band. Moreover, we marginally detected a phase offset of 33° ± 14° to the west of the planet. Both of these effects can be explained if TOI-1347 b is retaining a high-mean-molecular-weight atmosphere. Silicate clouds in this atmosphere could produce a high albedo, while partial cloud coverage may produce the observed phase offset to the West (see, e.g., Kepler-7 b, Demory et al. 2013). However, the data in hand are insufficient to definitively confirm the presence of an atmosphere on TOI-1347 b. Higher S/N follow-up observations with JWST are required.

We chose to adopt the same two-component simple harmonic oscillator (SHO) model for our RV model that we used to model stellar variability in photometry. We used the exoplanet (Foreman-Mackey et al. 2021) package to construct a Gaussian process (GP) stellar activity model with a kernel defined by the built-in RotationTerm parameterization. This kernel is a mixture of two SHOs analogous to the first and second cosine terms in Equation (1). We also tried a GP with the quasiperiodic kernel implemented in radvel (Fulton et al. 2018), both untrained and trained on the S-Indices. We found that the quasiperiodic kernel struggled to identify a single primary period, often jumping between 8 and 16 days. In order for the MCMC to converge, a reasonably strong prior (±0.5 day or less) had to be placed on one of these period solutions. We found this undesirable compared to the SHO kernel in exoplanet, which was able to lock onto the 16 day period even with wide priors of >10 days. We also found the final fit to be statistically indistinguishable whether the GP was trained on the S-Indices or not, so for our final model we opted for an untrained SHO GP.

For the full mathematical definition of the SHO GP kernel, the interested reader is directed to Foreman-Mackey et al. (2017). In brief, the kernel is parameterized by the GP amplitude (A_GP), the primary period of variability (P_GP), the quality factor of the primary mode Q₀, the difference in quality factors between the period and half-period modes (ΔQ₀), and lastly the fractional amplitude between the two modes (f). We followed the guidance of the exoplanet tutorials to choose the appropriate priors on these parameters, which are tabulated in Table 2.

Table 2. Transit and RV Parameters of the TOI-1347 System

Parameter	Symbol	Prior	Posterior Distribution	Unit
TOI-1347 b
Planet/Star Radius Ratio	Rp /R⋆	Uniform( − 10, 10) on $\mathrm{ln}({R}_{p}/{R}_{\star })$	0.02039 ± 0.00072	⋯
Time of Conjunction	Tc	Uniform(1682, 1683)	1682.71214 ± 0.00060	(BJD-2457000)
Orbital Period	Porb	Uniform( − 10, 10) on $\mathrm{ln}({P}_{\mathrm{orb}})$	0.84742346 ± 0.00000061	days
Orbital Inclination	iorb	Uniform( − 1, 1) on $\cos {i}_{\mathrm{orb}}$	79.5 ± 0.7	deg
Orbital Eccentricity	e	⋯	0 (fixed)	⋯
Impact Parameter	b	Derived	0.82 ± 0.03	⋯
Scaled Semimajor Axis	a/R⋆	Derived	4.43 ± 0.21	⋯
RV Semiamplitude	K	Normal(0, 50)	${7.74}_{-0.79}^{+0.80}$	m s−1
Planetary Radius	Rp	Derived	1.8 ± 0.1	R⊕
Planetary Mass	Mp	Derived	11.1 ± 1.2	M⊕
Bulk Density	ρ	Derived	${9.9}_{-1.7}^{+2.1}$	g cm−3
Equilibrium Temperature	Teq	Derived (AB = 0.7)	1400 ± 40	K
TOI-1347 c
Planet/Star Radius Ratio	Rp /R⋆	Uniform( − 10, 10) on $\mathrm{ln}({R}_{p}/{R}_{\star })$	0.0179 ± 0.0010	⋯
Time of Conjunction	Tc	Uniform(1678, 1679)	1678.5059 ± 0.0021	(BJD-2457000)
Orbital Period	Porb	Uniform( − 10, 10) on $\mathrm{ln}({P}_{\mathrm{orb}})$	4.841962 ± 0.000012	days
Orbital Inclination	i	Uniform( − 1, 1) on $\cos {i}_{\mathrm{orb}}$	87.5 ± 0.4	deg
Orbital Eccentricity	e	⋯	0 (fixed)	⋯
Impact Parameter	b	Derived	0.73 ± 0.08	⋯
Scaled Semimajor Axis	a/R⋆	Derived	14.18 ± 0.49	⋯
RV Semiamplitude	K	Normal(0, 50)	${1.08}_{-0.92}^{+0.91}\,(\lt 2.59\,\mathrm{at}\ 95 \% )$	m s−1
Planetary Radius	Rp	Derived	1.6 ± 0.1	R⊕
Planetary Mass	Mp	Derived	2.8 ± 2.3 ( < 6.4 at 95%)	M⊕
Bulk Density	ρ	Derived	${3.6}_{-3.0}^{+3.3}\,(\lt 4.1\,\mathrm{at}\ 95 \% )$	g cm−3
Equilibrium Temperature	Teq	Derived (AB = 0.1)	1000 ± 25	K
Stellar Parameters
Stellar Density	ρ⋆	Normal(2.25, 0.26)	2.30 ± 0.24	g cm−3
Limb Darkening	q1	Uniform(0, 1)	0.28 ± 0.20	⋯
Limb Darkening	q2	Uniform(0, 1)	0.35 ± 0.28	⋯
HIRES RV Jitter	${\sigma }_{\mathrm{HIRES}}$	HalfNormal(0, 10)	${4.05}_{-0.67}^{+0.73}$	m s−1
HARPS-N RV Jitter	σHARPS	HalfNormal(0, 10)	${8.20}_{-1.98}^{+2.67}$	m s−1
GP Amplitude	AGP	InverseGamma(1, 5)	${9.70}_{-1.06}^{+1.16}$	m s−1
Rotation Period	PGP	$\mathrm{Normal}(\mathrm{ln}(16),0.5)$ on $\mathrm{ln}({P}_{\mathrm{GP}})$	${16.16}_{-0.35}^{+0.35}$	days
Oscillator Quality Factor	Q0	HalfNormal(0, 2) on $\mathrm{ln}({Q}_{0})$	${1.49}_{-0.91}^{+0.74}$	⋯
Quality Factor Difference	ΔQ	Normal(0, 2) on $\mathrm{ln}({\rm{\Delta }}Q)$	${1.71}_{-2.33}^{+1.32}$	⋯
Fractional Amplitude	f	Uniform(0.1, 1)	${0.73}_{-0.24}^{+0.19}$	⋯
Instrumental RV Parameters
RV Offset (HIRES)	${\gamma }_{\mathrm{HIRES}}$	Normal(0, 200)	$-{2.12}_{-0.87}^{+0.90}$	m s−1
RV Offset (HARPS-N)	γHARPS	Normal(0, 200)	${6.29}_{-3.41}^{+3.35}$	m s−1

Download table as:  ASCIITypeset image

In practice, the Keplerian models are only parameterized by the RV semiamplitude K, which we place a wide, uninformative Gaussian prior on. We do not restrict K to positive values only, as such a prior can bias mass estimates to higher values, especially in the case of nondetections (Weiss & Marcy 2014). We imposed Gaussian priors on the periods and times of conjunction using the best-fit mean and standard deviation derived from our transit fits (Table 2), effectively fixing them to the tight transit constraints. Like the transit model, we fixed the eccentricity (e) and argument of periastron (ω) to zero. We included separate jitter terms for HIRES (${\sigma }_{\mathrm{HIRES}}$) and HARPS-N (σ_HARPS) and likewise separate RV offsets (${\gamma }_{\mathrm{HIRES}}$) and HARPS-N (γ_HARPS).

We initialized our model at the best-fit values from photometry, where applicable, and determined initial guesses for the RV semiamplitudes using exoplanet.estimate_semi_amplitude. We then solved for the MAP solution using scipy.optimize.minimize. The MAP parameters were then used as a seed for an MCMC exploration of the posterior. We employed a Hamiltonian Monte Carlo (HMC; Duane et al. 1987; Neal 2011) implemented in PyMC3 (Salvatier et al. 2016), specifically the No-U-Turn Sampler (NUTS; Hoffman & Gelman 2014). HMC and NUTS are generally more efficient than the traditional Metropolis–Hastings algorithm (Metropolis et al. 1953; Hastings 1970), resulting in well-mixed MCMC chains in far fewer samples. We sampled with NUTS in four parallel chains for 5000 "tuning" steps, which are discarded. The number of tuning steps was chosen to obtain an acceptance fraction near the target of 0.9, which balances the number of retained samples and efficient exploration of the posterior space. Each chain then collects 5000 samples for a total of 20,000 posterior samples. We ensured adequate statistical independence among these samples by requiring that the per-parameter $\hat{R}$ statistic (Vehtari et al. 2021) be <1.001.

Our best-fitting RV model is shown in Figure 7. The GP robustly recovers a primary period of 16.2 ± 0.3 days, even with a wide prior, independently verifying our assessment of the stellar rotation period from photometry. The USP TOI-1347 b is also robustly detected at consistent semiamplitudes, regardless of the activity model used (trained/untrained, quasiperiodic/SHO). The resulting mass of TOI-1347 b is 11.1 ± 1.2 M_⊕. The second planet, TOI-1347 c, is not detected in the RVs. We adopt an upper limit of <6.4 M_⊕ at 95% confidence. The residual rms to the combined Keplerian + GP model is 4.9 m s⁻¹ for HIRES and 7.5 m s⁻¹ for HARPS-N, which is (expectedly) similar to the fitted stellar jitter values (4.0 ± 0.7 m s⁻¹ and ${8.2}_{-2.0}^{+2.7}$ m s⁻¹ respectively), but larger than the per-measurement uncertainty of each instrument (2.5 and 1.9 m s⁻¹).

Zoom In Zoom Out Reset image size

Of note is the large jitter for the HARPS-N RVs. We suspect this is due to the GP being primarily conditioned on the HIRES RVs, resulting in poor predictive accuracy for activity during the timespan of the HARPS-N data, which occur about 1–2 rotation cycles before the HIRES data. As a result, the jitter term for HARPS-N is inflated to compensate. We tried to separate GPs for both data sets, sharing all hyperparameters except for the GP amplitude (A_GP). We found a nearly identical result (in fact, the best-fit HARPS-N jitter was higher) with statistically indistinguishable planet parameters, likely due to the relatively few HARPS-N data points. As a result, we adopt the single GP model but encourage further investigation into the nature of stellar activity on TOI-1347. As it stands, there are too few HARPS-N data points to condition independent GPs, but the single GP model is potentially overfitted to the HIRES points, reducing its out-of-sample predictive accuracy (Blunt et al. 2023).

TOI-1347 b is the largest (in both mass and radius) of the super-Earth USPs to date. It seems to be rocky in composition and similar in iron core mass fraction to the Earth. Figure 8 shows the two planets in the context of known exoplanets on the mass–radius diagram. ³² We modeled the core composition of TOI-1347 b assuming a simple two-layer model with an iron core and a silicate ("rock") mantle (Zeng et al. 2016). Our mass and radius measurements suggest an iron core mass fraction of 41% ± 27%, not far from Earth's 33% core mass fraction. TOI-1347 b joins a group of well-characterized USPs (Dai et al. 2019) that are consistent with an Earth-like composition.

Zoom In Zoom Out Reset image size

With a core mass larger than 10 M_⊕, TOI-1347 b, along with the USPs TOI-1075 b (Essack et al. 2023) and HD 20329 b (Murgas et al. 2022), are close to the theoretical limit for runaway accretion (see, e.g., Rafikov 2006). How did these planets evade runaway accretion and not become gas giants? Lee (2019) and Chachan et al. (2021) both noted that the local hydrodynamic conditions, the envelope opacity, and the timescale of core assembly relative to disk dissipation could all contribute to quenching runaway accretion. As more of these systems are discovered, population-level analyses may shed light on which planets are able to evade runaway and which grow into gas giants.

TOI-1347 b also pushes the efficacy of photoevaporation to its limit. At >10 M_⊕, the outflowing atmosphere has to overcome a deep gravitational potential well. On the other hand, the temperature of atmospheric outflows is likely capped below 10⁴ K due to strong radiative cooling at higher temperatures (see, e.g., Murray-Clay et al. 2009). An order-of-magnitude comparison reveals that the thermal sound speed (∼10 km s⁻¹ at 10⁴ K) may not overcome the escape velocity of the planet (∼25 km s⁻¹ at 10 M_⊕ and 1.9 R_⊕), preventing bulk hydrodynamic outflow (i.e., photoevaporation). Previous models showed that photoevaporation is significantly quenched on planets with heavier cores (≳6 M_⊕; see, e.g., Owen & Wu 2017; Wang & Dai 2018). Planets like TOI-1347 b are therefore important test cases to understand the limit of both photoevaporation and core-powered mass loss (Ginzburg et al. 2018; Gupta & Schlichting 2019). Our mass and radius measurements disfavor the presence of a thick H/He envelope (Figure 8). Did TOI-1347 b have, and then lose, a primordial H/He envelope? Or, could TOI-1347 b have formed near to its present-day scorching orbit without ever acquiring a substantial atmosphere? We encourage further investigation on this question.

Even though TOI-1347 b has a mass and radius that suggest an Earth-like bulk composition, we cannot rule out the presence of a heavy-mean-molecular-weight atmosphere. Lopez (2017) showed that the largest of the nongiant USPs (such as 55 Cnc e) can hold onto high-metallicity atmospheres even in the presence of strong stellar radiation. Such an atmosphere would only marginally inflate the planet's radius (the scale height of a CO₂-based atmosphere is about 11 km, while planetary radii are ∼10,000 km).

In fact, our tentative phase curve (3σ) and secondary eclipse (2σ) detections of TOI-1347 b point to a nonzero albedo and a possible phase offset to the west. These features could indicate the presence of an atmosphere at least partially covered by reflective silicate clouds. It may be the case that the deep gravitational wells of the most massive super-Earth USPs are sufficient to cling to such an atmosphere and resist atmospheric loss mechanisms. R. Hu et al. (2024, in preparation) show that their JWST NIRCam and MIRI observations of 55 Cnc e, a similar massive USP (0.73 day orbit; 9M_⊕), can only be explained if the planet still has a CO- or CO₂-based atmosphere. Future observations with JWST might uncover a similar story for TOI-1347 b (TSM =19, ESM =5.6, using the equations of Kempton et al. 2018).

Alternatively, the high-amplitude TESS phase curve of TOI-1347 b may be a consequence of outgassed Na emission on the hot dayside of the planet. Zieba et al. (2022) showed that this emission can explain the phase curve of the lava world K2-141 b, which has been observed in both the Kepler passband and Spitzer-4.5 μm passband (Malavolta et al. 2018). In their analysis, they found that the two phase curves are inconsistent with a blackbody model, with the visible-light phase curve having a higher amplitude than expected. Similarly high visible-light phase-curve amplitudes have been reported for the lava worlds 55 Cnc e (Kipping & Jansen 2020) and Kepler-10 b (Batalha et al. 2011; Rouan et al. 2011) although the latter has not yet been observed in the infrared, which would test blackbody emission. If Na emission is responsible for the observations of TOI-1347 b reported here, it may more easily explain the tentative phase-curve offset than reflective clouds, which would need to be nonuniformly distributed across the planet.